
*==============================================================================*
*      Used Vehicle Emissions Rate and Miles Driven, by Vintage and Age        *
*==============================================================================*



u "dataSTATA/combined/combined_smogcheck_colorado.dta", clear
keep if age >= 4

* age in 2-year categories
replace age = 4  if age == 5
replace age = 6  if age == 7
replace age = 8  if age == 9
replace age = 10 if age == 11
replace age = 12 if age == 13
replace age = 14 if age == 15
replace age = 16 if age == 17
replace age = 18 if age == 19
replace age = 20 if age == 21
replace age = 22 if age == 23
replace age = 24 if age == 25
replace age = 26 if age == 27
replace age = 28 if age == 29
replace age = 30 if age >= 31 & age < .


g       group = 1 if model_year >= 2004 & model_year <= 2009
replace group = 2 if model_year >= 1998 & model_year <= 2003
replace group = 3 if model_year >= 1992 & model_year <= 1997
replace group = 4 if model_year <= 1991 

* exclude miniscule cell sizes
keep if (group == 1 & age <= 10) | (group == 2 & age <= 16) | (group == 3 & age <= 22) | (group == 4 & age >= 6)

preserve
	collapse (mean) emissions_* annualVmt , by(group age) fast

	tw	(connected emissions_used_CO age if group == 4, msymb(C))	///
		(connected emissions_used_CO age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_CO age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_CO age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("CO Emissions (Grams per Mile)")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "Model years <1992" 2 "Model years 1992-1997" 3 "Model years 1998-2003" 4 "Model years 2004-2009")) ///
		ylab(1 2 4 8 16 32) ///
		xlab(4 8 12 16 20 24 30) ///
		saving("results/figures/f6_cohorts/CO.gph", replace)
	graph export "results/figures/f6_cohorts/CO.wmf", replace
	graph export "$overleaf/figures/f6_cohorts/CO.eps", replace
		
	tw	(connected emissions_used_HC age if group == 4, msymb(C))	///
		(connected emissions_used_HC age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_HC age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_HC age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("HC Emissions (Grams per Mile)")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "Model years <1992" 2 "Model years 1992-1997" 3 "Model years 1998-2003" 4 "Model years 2004-2009")) ///
		ylab(0.05 0.1 0.2 0.4 0.8 1.6 3.2) ///
		xlab(4 8 12 16 20 24 30) ///
		saving("results/figures/f6_cohorts/HC.gph", replace)
	graph export "results/figures/f6_cohorts/HC.wmf", replace
	graph export "$overleaf/figures/f6_cohorts/HC.eps", replace

	tw	(connected emissions_used_NOX age if group == 4, msymb(C))	///
		(connected emissions_used_NOX age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_NOX age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_NOX age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("NO{subscript:X} Emissions (Grams per Mile)")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "Model years <1992" 2 "Model years 1992-1997" 3 "Model years 1998-2003" 4 "Model years 2004-2009")) ///
		ylab(0.25 0.5 1 2) ///
		xlab(4 8 12 16 20 24 30) ///
		saving("results/figures/f6_cohorts/NOX.gph", replace)
	graph export "results/figures/f6_cohorts/NOX.wmf", replace
	graph export "$overleaf/figures/f6_cohorts/NOX.eps", replace

	tw	(connected emissions_used_CO2 age if group == 4, msymb(C))	///
		(connected emissions_used_CO2 age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_CO2 age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_CO2 age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("CO{subscript:2} Emissions (Grams per Mile)")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "Model years <1992" 2 "Model years 1992-1997" 3 "Model years 1998-2003" 4 "Model years 2004-2009")) ///
		ylab(200 400 800) ///
		xlab(4 8 12 16 20 24 30) ///
		saving("results/figures/f6_cohorts/CO2.gph", replace)
	graph export "results/figures/f6_cohorts/CO2.wmf", replace
	graph export "$overleaf/figures/f6_cohorts/CO2.eps", replace

	tw	(connected annualVmt age if group == 4, msymb(C))	///
		(connected annualVmt age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected annualVmt age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected annualVmt age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("Annual VMT ('000s)")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "Model years <1992" 2 "Model years 1992-1997" 3 "Model years 1998-2003" 4 "Model years 2004-2009")) ///
		ylab(20 15 10 7.5 5) ///
		xlab(4 8 12 16 20 24 30) ///
		saving("results/figures/f6_cohorts/VMT.gph", replace)
	graph export "results/figures/f6_cohorts/VMT.wmf", replace
	graph export "$overleaf/figures/f6_cohorts/VMT.eps", replace

restore


***** version with figure subtitles for slides

preserve
	collapse (mean) emissions_* annualVmt , by(group age)

	tw	(connected emissions_used_CO age if group == 4, msymb(C))	///
		(connected emissions_used_CO age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_CO age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_CO age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("Log Emissions")		 ///
		xtit("Age (Years)")				 ///
		legend(off) ///
		ylab(1 2 4 8 16 32) ///
		xlab(4 8 12 16 20 24 30) ///
		subtit("CO") ///
		saving("results/figures/f6_cohorts/CO_forslides.gph", replace)
		
	tw	(connected emissions_used_HC age if group == 4, msymb(C))	///
		(connected emissions_used_HC age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_HC age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_HC age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("Log Emissions")		 ///
		xtit("Age (Years)")				 ///
		legend(off) ///
		ylab(0.05 0.1 0.2 0.4 0.8 1.6 3.2) ///
		xlab(4 8 12 16 20 24 30) ///
		subtit("HC") ///
		saving("results/figures/f6_cohorts/HC_forslides.gph", replace)

	tw	(connected emissions_used_NOX age if group == 4, msymb(C))	///
		(connected emissions_used_NOX age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_NOX age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_NOX age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("Log Emissions")		 ///
		xtit("Age (Years)")				 ///
		legend(off) ///
		ylab(0.25 0.5 1 2) ///
		xlab(4 8 12 16 20 24 30) ///
		subtit("NOX") ///
		saving("results/figures/f6_cohorts/NOX_forslides.gph", replace)

	tw	(connected emissions_used_CO2 age if group == 4, msymb(C))	///
		(connected emissions_used_CO2 age if group == 3, msymb(S) lpattern(dash))	 ///
		(connected emissions_used_CO2 age if group == 2, msymb(Dh) lpattern(shortdash))	 ///
		(connected emissions_used_CO2 age if group == 1, msymb(Th) lpattern(longdash)),	 ///
		graphr(color(white))				 ///
		yscale(noline log)				 ///
		ytit("Log Emissions")		 ///
		xtit("Age (Years)")				 ///
		legend(order(1 "MY <1992" 2 "MY 92-97" 3 "MY 98-03" 4 "MY 04-09")) ///
		ylab(200 800) ///
		xlab(4 8 12 16 20 24 30) ///
		subtit("CO2") ///
		saving("results/figures/f6_cohorts/CO2_forslides.gph", replace)

restore

graph combine "results/figures/f6_cohorts/CO_forslides.gph" ///
	"results/figures/f6_cohorts/HC_forslides.gph" ///
	"results/figures/f6_cohorts/NOX_forslides.gph" ///
	"results/figures/f6_cohorts/CO2_forslides.gph", ///
	graphr(color(white)) 
graph export "results/figures/f6_cohorts/COHCNOXCO2_forslides.pdf", replace
	
